---
title: "IMCE and VIMP estimation"
author: "Dai Yamao"
date: "`r format(Sys.time(), '%Y-%m-%d')`"
output:
  html_document: default
  pdf_document: default
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, fig.width = 10, fig.height = 5)
rm(list=ls())

require(tidyverse)
require(cjoint)
#library(cregg)
library(scales)
library(patchwork)
library(extrafont)
require(memisc) # write_html()
library(cjbart)
library(Cairo)

dat <- read_csv("data/cj_data.csv")
dat <- na.exclude(dat)

dat <- dat |>
  dplyr::select(selected, General_description, Internal_external,
                Religious_customary, Support_criticism, respondentIndex,
                Q2_1, Q2_3,
                Q14,
                Gender, D8, D6, D7,
                Q4, Q5, Q6, Q7,
                Q4_1, Q4_2, Q4_3)
str(dat)
dat[,c(2:5)] <- data.frame(lapply(dat[,c(2:5)],as.factor))

# covariate for VIMP
dat <- dat |> 
  mutate(islamism = Q2_1,
         tribalism = Q2_3,
         Patronege_rel_daily = Q4,
         Patronege_rel_job = Q5,
         Patronege_trib_daily = Q6,
         Patronege_trib_job = Q7,
         SLA = Q4_1,
         Sadr = Q4_2,
         Fatah = Q4_3)


# Islamism
dat <- dat |> 
  mutate(Islamism = case_when(Q2_1 >= 67 ~ "High",
                              Q2_1 <= 66 & Q2_1 >= 34 ~ "Middle",
                              Q2_1 <= 33 ~ "Low"))
dat$Islamism <- as.factor(dat$Islamism)

# Tribalism
dat <- dat |> 
  mutate(Tribalism = case_when(Q2_3 >= 67 ~ "High",
                               Q2_3 <= 66 & Q2_3 >= 34 ~ "Middle",
                               Q2_3 <= 33 ~ "Low"))
dat$Tribalism <- as.factor(dat$Tribalism)

# Child marriage
dat <- dat |> 
  mutate(Childe_marriage = case_when(Q14 >= 3 ~ "Agree",
                                     Q14 == 2 ~ "Neither",
                                     Q14 == 1 ~ "Disagree"))
dat$Childe_marriage <- as.factor(dat$Childe_marriage)


# gender
dat <- dat |> 
  mutate(Gender = case_when(Gender == 1 ~ "Male",
                            Gender == 2 ~ "Female"))
dat$Gender <- as.factor(dat$Gender)

# urban/rural
dat <- dat |> 
  mutate(Urban_rural = case_when(D8 == 1 ~ "Urban",
                                 D8 == 2 ~ "Rural"))
dat$Urban_rural <- as.factor(dat$Urban_rural)

# sect
dat <- dat |> 
  mutate(Sect = case_when(D6 == 1 ~ "Sunni",
                          D6 == 2 ~ "Shia",
                          D6 == 3 ~ "Christian",
                          D6 == 4 ~ "Kurd",
                          D6 == 5 ~ "others"))
dat$Sect <- as.factor(dat$Sect)

# patronage_religious_daily_necessities
dat <- dat |> 
  mutate(ptr_rel_daily = case_when(Q4 >= 4 ~ "High",
                                   Q4 <= 2 ~ "Low",
                                   Q4 == 3 ~ "other"))
dat$ptr_rel_daily <- as.factor(dat$ptr_rel_daily)

# patronage_religious_job
dat <- dat |> 
  mutate(ptr_rel_job = case_when(Q5 >= 4 ~ "High",
                                 Q5 <= 2 ~ "Low",
                                 Q5 == 3 ~ "other"))
dat$ptr_rel_job <- as.factor(dat$ptr_rel_job)

# patronage_Tribe_daily_necessities
dat <- dat |> 
  mutate(ptr_trb_daily = case_when(Q6 >= 4 ~ "High",
                                   Q6 <= 2 ~ "Low",
                                   Q6 == 3 ~ "other"))
dat$ptr_trb_daily <- as.factor(dat$ptr_trb_daily)

# patronage_Tribe_job
dat <- dat |> 
  mutate(ptr_trb_job = case_when(Q7 >= 4 ~ "High",
                                 Q7 <= 2 ~ "Low",
                                 Q7 == 3 ~ "other"))
dat$ptr_trb_job <- as.factor(dat$ptr_trb_job)


# Shia parties supporter
dat <- dat |> 
  mutate(Shia_party_supporter = (Q4_1 + Q4_2 + Q4_3)/3) |>
  mutate(Shia_party_supporter = ifelse(Shia_party_supporter >= 7, "Shia Party Supporter", "others"))
dat$Shia_party_supporter <- as.factor(dat$Shia_party_supporter)


# region south
dat <- dat |> 
  mutate(South = ifelse(D7 >= 12, "South", "other"))
dat$South <- as.factor(dat$South)

# rural south
dat <- dat |>
  mutate(rural_south = ifelse(South == "South" & Urban_rural == "Rural", "Rural_south", "Other"))
dat$rural_south <- as.factor(dat$rural_south)


# region west
dat <- dat |> 
  mutate(West = ifelse(D7 >= 4 & D7 <= 8, "West", "other"))
dat$West <- as.factor(dat$West)


# South Islamism
dat <- dat |>
  mutate(South_Islamism = ifelse(Islamism == "High" & South == "South", "South_Islamism", "Other"))
dat$South_Islamism <- as.factor(dat$South_Islamism)

# West Tribalism
dat <- dat |>
  mutate(West_tribalism = ifelse(Tribalism == "High" & West == "West", "West_Tribalism", "Other"))
dat$West_tribalism <- as.factor(dat$West_tribalism)

# South Tribalism
dat <- dat |>
  mutate(South_tribalism = ifelse(Tribalism == "High" & South == "South", "South_Tribalism", "Other"))
dat$South_tribalism <- as.factor(dat$South_tribalism)

# rural Islamism
dat <- dat |>
  mutate(rural_Islamism = ifelse(Islamism == "High" & Urban_rural == "Rural", "Rural_Islamism", "Other"))
dat$rural_Islamism <- as.factor(dat$rural_Islamism)

# rural Tribalism
dat <- dat |>
  mutate(rural_tribalism = ifelse(Tribalism == "High" & Urban_rural == "Rural", "Rural_tribalism", "Other"))
dat$rural_tribalism <- as.factor(dat$rural_tribalism)

```


# IMCE
```{r, fig.height = 8, dpi = 300, warning = FALSE}
dat <- dat |>
  mutate(Support_criticism = dplyr::recode(Support_criticism,
                      "Legal fragmentaion" = "Legal_fragmentaion",
                      "Protect poor women" = "Protect_poor_women",
                      "Social order" = "Social_order",
                      "Violation to women's right" = "Violate_Women_right"))

set.seed(89)

# modeling
cj_model <- cjbart(data = dat,
                   Y = "selected",
                   type = "rating",
                   id = "respondentIndex")

# IMCE estimation
het_effect <- IMCE(data = dat,
                   model = cj_model,
                   attribs = c("General_description", "Religious_customary", "Support_criticism","Internal_external"),
                   ref_levels = c("De fact legalization", "Religious laws", "Legal fragmentaion","Domestic criticisms"),
                   cores = 32,
                   method = "bayes",
                   keep_omce = FALSE)
summary(het_effect)
plot(het_effect)

plot(het_effect,
     covar = c("islamism"),
     plot_levels = c("Right to choose", "Tribal customs","Violate_Women_right"))

plot(het_effect,
     covar = c("tribalism"),
     plot_levels = c("Right to choose", "Tribal customs","Violate_Women_right"))

# plot
plot(het_effect,
     covar = c("islamism"))

plot(het_effect,
     covar = c("islamism"))

plot(het_effect,
     covar = c("tribalism"))

plot(het_effect,
     covar = c("SLA"))

plot(het_effect,
     covar = c("Sadr"))

plot(het_effect,
     covar = c("Fatah"))

plot(het_effect,
     covar = c("Patronege_rel_daily"))

plot(het_effect,
     covar = c("Patronege_trib_daily"))

plot(het_effect,
     covar = c("Urban_rural"))
     
```

# IMCE plot
```{r, fig.height = 5, dpi = 300, warning = FALSE}
H1a <- plot(het_effect,
     covar = c("islamism"),
     plot_levels = c("Violate_Women_right", "Protect_poor_women"))
H1a
ggsave(filename = "result/H1a.pdf", 
       plot = H1a, 
       dpi = 300, 
       width = 7, 
       height = 5)

H1b <- plot(het_effect,
     covar = c("Shia_party_supporter"),
     plot_levels = c("International criticism", "Domestic criticism"))
H1b
ggsave(filename = "result/H1b.pdf", 
       plot = H1b, 
       dpi = 300, 
       width = 7, 
       height = 5)

H1c <- plot(het_effect,
     covar = c("Patronege_rel_daily"),
     plot_levels = c("Right to choose", "Tribal customs", "Evaluation"))
H1c
ggsave(filename = "result/H1c.pdf", 
       plot = H1c, 
       dpi = 300, 
       width = 7, 
       height = 5)

H2a <- plot(het_effect,
     covar = c("tribalism"),
     plot_levels = c("Social_order", "Evaluation"))
H2a
ggsave(filename = "result/H2a.pdf", 
       plot = H2a, 
       dpi = 300, 
       width = 7, 
       height = 5)

H2b <- plot(het_effect,
     covar = c("Patronege_trib_daily"),
     plot_levels = c("Right to choose", "Violate_Women_right", "Social_order", "Protect_poor_women"))
H2b
ggsave(filename = "result/H2b.pdf", 
       plot = H2b, 
       dpi = 300, 
       width = 7, 
       height = 5)
```

# VIMP
```{r, fig.height = 7, fig.width = 8, dpi = 300, warning = FALSE}
# random forest model
vimp <- het_vimp(het_effect,
                 covars = c("islamism", "tribalism", 
                            "SLA",
                            "Sadr",
                            "Fatah",
                            "Patronege_rel_daily",
                            "Patronege_trib_daily", 
                            "Shia_party_supporter"),
                 cores = 32)

plot(vimp,
     covars = c("islamism", "tribalism", 
                "SLA",
                "Sadr",
                "Fatah",
                "Patronege_rel_daily",
                "Patronege_trib_daily", 
                "Shia_party_supporter"))

# plot heatmap
covar_lookup <- data.frame(covar_name = c("islamism", "tribalism", 
                                          "SLA",
                                          "Sadr",
                                          "Fatah",
                                          "Patronege_rel_daily",
                                          "Patronege_trib_daily", 
                                          "Shia_party_supporter"),
                           covar_label = c("Islamism", "Tribalism", 
                                           "SLA",
                                           "Sadr",
                                           "Fatah",
                                           "Religious Patronege",
                                           "Tribal Patronege", 
                                           "Shia_party_supporter"))

final_results <- vimp$results %>% 
  left_join(covar_lookup, by = c("covar" = "covar_name")) %>% 
  mutate(lab_text = ifelse(importance < 10,"",round(importance,2)),
         Attribute = case_when(Attribute == "General_description" ~ "Dscrp",
                               Attribute == "Religious_customary" ~ "Rel/trib",
                               Attribute == "Support_criticism" ~ "Evaluation",
                               Attribute == "Internal_external" ~ "Internal or external"))

p20 <- ggplot(final_results, aes(x = covar_label, y = Level, fill = importance)) +
  facet_grid(Attribute ~ ., space = "free", scales = "free", switch = "y") +
  geom_tile() +
  scale_fill_gradient(low="white", high="firebrick2") +
  labs(x = "Subject covariate", y = "Attribute-level", fill = "Importance") +
  theme(text = element_text(size = 14),
        axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1))
p20

ggsave(filename = "result/VIMP.pdf", 
       plot = p20, 
       dpi = 300, 
       width = 8, 
       height = 7)

```
